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Reverse Monte Carlo (RMC) is an algorithm that incorporates stochastic modification of the ac¬ 
tion as part of the process that updates the fields in a Monte Carlo simulation. Such update moves 
have the potential of lowering or eliminating potential barriers that may cause inefficiencies in ex¬ 
ploring the field configuration space. The highly successful Cluster algorithms for spin systems 
can be derived from the RMC framework. In this work we apply RMC ideas to pure gauge theory, 
aiming to alleviate the critical slowing down observed in the topological charge evolution as well 
as other long distance observables. We present various formulations of the basic idea and report 
on our numerical experiments with these algorithms. 
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1. Introduction 

Critical slowing down is one of the main hindrances to reliable computation with Monte Carlo 
simulations in Lattice Field Theory. This is the phenomenon of growing computational costs when 
approaching a theory’s critical point because of large autocorrelation times. The severity of this 
slow down is dependent on the algorithm used and the observable of interest. In order to obtain 
reliable statistics, all Monte Carlo simulations must have a much longer simulation time than the 
autocorrelation time of the measured observable. 

In Lattice QCD, when approaching the continuum limit, the topological charge develops a very 
large autocorrelation time ||T]]. This is because of the large barriers between topological sectors that 
develop at small lattice spacing, thereby reducing the tunneling rate between these sectors. Large 
autocorrelation times for the topological charge as well as for other long distance observables, 
represent a major challenge for current QCD simulations. 

Recently, open boundary conditions in the Euclidean time direction instead of periodic have 
been proposed 0, as a solution to this problem. This allows topological charge to flow through 
the boundaries of the lattice and alleviate critical slowing down. In recent work the effectiveness 
of open boundary conditions has been studied extensively with Hybrid Monte Carlo (HMC) simu¬ 
lations and flow of topological charge has been shown to be a diffusive process [^]. 

Although open boundary conditions decrease the autocorrelation time at fine lattice spacing, 
translation invariance in the time direction is lost in these simulations. In certain cases it is prefer¬ 
able to use periodic boundaries. This gives motivation for new algorithms that decrease autocorre¬ 
lation time yet retain periodic boundary conditions to be explored. In this work we explore possible 
modifications to the Monte Carlo update process that may reduce autocorrelation times. One such 
algorithm is Reverse Monte Carlo (RMC) 0 that was found to be very effective in scalar field 
theories in two dimensions. 

2. Reverse Monte Carlo 

Reverse Monte Carlo (RMC) is an algorithm that apart from updating the fields of a theory, 
also introduces local stochastic modification of the action itself. We write the action for the system 
we are simulating as 

S = £^(n), (2.1) 

n 

where n is an index of the lattice site and .S^{n) is the local action density. We introduce local 
modifications of the action density .ifi(n) that can occur with probability P{n) on each site n given 
by 

P{n) = = - - -, (2.2) 

where C is a normalization constant that ensures the switch probability is P{n) < 1 for all possible 
field values on site n. One such convenient choice is 

C = , (2.3) 

where the maximum is taken over all possible values of the fields in the action. Note that the 
choice of C is not unique and any choice that ensures P{n) < 1 is sufficient. On a given field 
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configuration during the simulation, local changes of the action are proposed and accepted with 
probability P{n). However, in order to have a Monte Carlo update that converges to the desired 
probability distribution dictated by the original action, the action density of sites that do not switch 
needs to also be modified as shown in [0]. The resulting action density after switching becomes 


^'{n) 


££s (n) sites that switched 

^{n) — ln(l — P{n)) sites that did not switch 


(2.4) 


Subsequently the fields are updated, but now using the new modified action for the system. Any 
update method of the fields can be incorporated, however maximal efficiency can be achieved if 
specific features of the modified action are exploited as was done in 0. For details on the proof 
of validity of this method the reader is referred to 0. It turns out that this approach represents 
a generalization to cluster algorithms for spin systems and in fact the original Swendsen-Wang 
cluster algorithm |^] for the Ising model can be derived in this framework. 


3. Application of RMC to QCD 


The application of RMC to Lattice QCD is straight forward. Our example is with pure gauge 
theory and the Wilson gauge action. Incorporation of fermions is straight forward given that our 
basic field update algorithm is Hybrid Monte Carlo (HMC) and thus it is omitted here. Furthermore, 
the issue of large autocorrelations of the topological charge arises in pure gauge theory as well. 

The Wilson gauge action is written as 

Sg [U]=P £ ^ReTr [l - U^y{n)] = £ A^v(«), (3.1) 

n,g,v 


where U^v{n) is the plaquette and A^y{n) is the action density on site n. The pure gauge partition 
function is 



(3.2) 


The local modification to the action density we are proposing is to introduce a plaquette and 
site dependent coupling The coupling constant can be allowed to fluctuate between the 

original value jS and a lower value a. If an appreciable fraction of plaquettes acquire a small 
coupling constant a, one can hope that updates of the gauge fields with this modified action can 
allow topological charge fluctuations that are suppressed with the original action. 

The probability for a plaquette to switch is dependent on the difference of unswitched and 
switched terms 

^ ^ (3 3 ) 


with 


_ ^max[{ji-a)A^y{n)] 


(3.4) 


Due to SU(3) being a compact group, the trace of the plaquette is bounded, hence there exists a 
unique value for C. 

The new action to be used in the subsequent update is 


[^] ~ ^ Pgv{tT-) A^v{n) ^ ln(l Pp^(j{n)) 

n,g,v neA^,p.a 


(3.5) 
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where JV is the set of sites which did not switch, and retained the original coupling j3. In our 
application we use HMC for the subsequent field update, hence the inclusion of fermions is straight 
forward. 


3.1 Continuous Beta 

Direct application of RMC to QCD results in small switch rates, especially if a large change 
in P is allowed. For this reason one may want to modify the algorithm in a way that gives a range 
of P values. Following this idea, we introduce a new method that allows for a continuous variation 
of the coupling constant in the interval [oc,p]. In this case the switch probability is given by: 


resulting in a new action with a new field iSav to be simulated given by 




U,Pjiv — ^ In 


1 




,(/3—o:)Af,v(f2) _ ^ 


(3.7) 


We call this algorithm “Continuous Beta”. Unlike traditional RMC, Continuous Beta gives a log 
term for every plaquette. The advantage of this method is that a can be set to 0 and plaquettes are 
more likely to switch close to this lower bound. The logarithmic piece that accompanies each in¬ 
teraction must be accounted for however, and in general introduces nontrivial tunneling properties. 
Note that in this case the chosen probability distribution for the coupling constant integrates to one 
i.e. 

/•/3 , 

/ p{Pgv{n)) dP^y{n) = 1. (3.8) 

J a 

There are other possibilities for implementing a continuous coupling algorithm, which are 
currently under investigation. The validity of this algorithm will be explained in the next section 
where an alternative proof for the RMC algorithm is also presented. 


4. Path Integral Interpretation 

The authors in 0 proved that Reverse Monte Carlo respects detailed balance of the original 
theory. When applying RMC to Lattice Field Theory, the same conclusion can be reached using a 
simple argument based on the path integral that defines the correlation functions to be simulated. 
The fluctuating coupling constants can be thought of as a new site dependent field. The action 
of this new field has been construcfed in such way that after integrating over this new field fhe 
original pafh integral with fixed coupling consfant is recovered. Therefore, simulating the new 
augmented theory results in the same correlation functions as the original theory. To be more 
precise, using Eq. ^ in the case of the Continuous Beta algorithm, the following identity holds for 
the partition function 

[ D[U] = [ D[U] (4.1) 

J J Ja yi 
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which results in 



(4.2) 


Hence, correlation functions of the original theory can be written as 


i^) = ^I D[U]D[p^y] ff{U) 


(4.3) 


where (^([/) is the desired correlation function of gauge fields. From these identities it is clear that 
the proposed update algorithm is nothing but a normal update algorithm for the new augmented 
theory that first updates the coupling constant field and subsequently updates the gauge fields. 

For RMC fhe same formalism holds with the distinction that the integral over the fluctuating 
coupling constant has to be replaced with a sum over two distinct values of j3. In this case the sum 
of probabilities for having the coupling constant to be j8 and probability of switching to a add up 
to one. Hence the manipulations done for the Continuous Beta case also hold. Furthermore, this 
point of view gives an easy way to generalize RMC to the case of multiple distinct values for the 
coupling constant. In addition, the path integral formalism makes it easy to design other update 
algorithms with fluctuating terms in the action. 

Now that we have established the correctness of this class of algorithms (RMC and the Con¬ 
tinuous Beta generalization), we may attempt to tune the update algorithm of the augmented ac¬ 
tion. One idea we tried was to introduce a freezing parameter into the update step of the coupling 
constants. Instead of switching the coupling on every update, we froze it for K steps. Setting a 
moderately high K ensures that the simulation has enough time to take advantage of potentially 
lowered barriers between topological sectors. From the work done in [^, it is observed that the 
topological charge fluctuations obey a simple diffusion and tunneling model. Our action modifica¬ 
tions may enhance the tunneling rates locally, however enough simulation time needs to be given 
for the diffusion process to propagate topological fluctuations. Careful study of the dependence of 
the algorithm on K is needed in order to construct an efficient algorithm. 

5. Simulation Results 

We performed numerical experiments on a 32'^ lattice with j3 = 6.179 with the pure gauge 
Wilson action. The simulation algorithm of choice was HMC with a fourth order integrator. The 
integrator parameters used resulted in an acceptance rate of « 90%. 

The Continuous Beta simulations utilized a lower bound for the coupling constant of a = 0 
and a = 5. The basic RMC algorithm was run with a modified coupling of a = 5. Additionally, 
separate experiments with RMC a = 5 that froze the stochastic modification for 100 or 200 steps 
were also conducted. We computed expectation values and autocorrelation times for the plaquette 
and smeared plaquette, and topological charge using the above simulation algorithms as well as the 
basic HMC algorithm. The pure HMC simulation is marked the "Base" case in subsequent tables. 
The smearing method used to compute topological charge and the smeared plaquette was gradient 
flow with a flow time of T = 2 in lattice units. 
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Table 1: Plaquette 10,000 Trajectories 


Algorithm 

Value 

Integrated Autocorrelation 

Base 

0.6116998(3) 

3.26(3) 

Continuous Beta 0 

0.611699(2) 

4.2(3) 

Continuous Beta 5 

0.611699(3) 

3.9(4) 

Reverse Monte Carlo 5 

0.611699(2) 

3.5(2) 

Reverse Monte Carlo 5 Freeze 100 

0.611711(4) 

5.1(5) 

Reverse Monte Carlo 5 Freeze 200 

0.611704(4) 

5.0(5) 


Table 2: Smeared Plaquette 10,000 Trajectories 


Algorithm 

Value 

Integrated Autocorrelation 

Base 

0.9989961(15) 

66(19) 

Continuous Beta 0 

0.9989934(8) 

43.(8) 

Continuous Beta 5 

0.9989935(14) 

58(16) 

Reverse Monte Carlo 5 

0.9989955(9) 

79(15) 

Reverse Monte Carlo 5 Freeze 100 

0.9989964(16) 

78(2) 

Reverse Monte Carlo 5 Freeze 200 

0.9989973(15) 

64(18) 


Table 3: Topological Charge 10,000 Trajectories 


Algorithm 

Value 

Integrated Autocorrelation 

Base 

1.02 ± 1.05 

320(140) 

Continuous Beta 0 

2.63 ± 1.06 

420(180) 

Continuous Beta 5 

1.80 ± 1.14 

260(110) 

Reverse Monte Carlo 5 

-0.23 ± 0.7 

190(80) 

Reverse Monte Carlo 5 Freeze 100 

-0.69 ± 1.05 

256(110) 

Reverse Monte Carlo 5 Freeze 200 

1.23 ± 1.54 

530(250) 


The analysis of the integrated autocorrelation time was carried out with the tools released by 
the alpha collaboration [Q]. Our results so far are largely inconclusive due to low statistics. We have 
taken a few of the algorithm variants to 100,000 trajectories, however the data is still not precise 
enough to give any definitive results. 

6. Conclusion 

Although tests so far have not indicated a dramatic improvement of topological freezing, the 
low statistics of our data prevents us from having firm conclusions. There is a vast parameter space 
unexplored and more extensive investigation on the dependence of autocorrelation times as func¬ 
tions of the switched coupling and the specifics of the simulation algorithm need to be made. The 
the path integral interpretation of these algorithms provides a simple method for generalizing these 
algorithms, possibly allowing us to find an optimal modification that improves the autocorrelation 
time of topological charge. We are currently investigating such generalizations that admit action 
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modifications at larger length scales, comparable to the longest correlation length in the theory. 
Such an approach may be more likely to produce a favorable outcome closer to the continuum 
limit. Additionally, we are also extending our statistics. An analysis with about a million trajecto¬ 
ries will be available soon. 
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